######################
##### PREPARATION ####
######################
# Clear environment
rm(list=setdiff(ls(), "path_to"))

# Load data
s = readRDS(path_to("wide"))
s = s %>% filter(!is.na(code_effect_simple))


######################
#### DERIVE DATA #####
######################
# Main estimates
df = rbind(
  s %>% group_by(vig, code_effect_simple) %>%
    summarise(
      n = n(),
      wgt = sum(wgt),
      nospeeder = sum(Duration__in_seconds_ > quantile(s$Duration__in_seconds_, 0.2)),
      consumer = sum(consumer_vig),
      strictconseq = sum(wtp_ratio_class == 0, na.rm=T),
      weakconseq = sum(wtp_ratio < 1, na.rm=T)
    ) %>%
    mutate(
      freq_main = n / sum(n, na.rm=T),
      freq_wgt = wgt / sum(wgt, na.rm=T),
      freq_nospeeders = nospeeder / sum(nospeeder),
      freq_consumer = consumer / sum(consumer),
      freq_strictconseq = strictconseq / sum(strictconseq),
      freq_weakconseq = weakconseq / sum(weakconseq)
  ) %>% filter(code_effect_simple != "dampening"),
  s %>% group_by(vig, code_effect_simple=code_effect) %>%
    summarise(
      n = n(),
      wgt = sum(wgt),
      nospeeder = sum(Duration__in_seconds_ > quantile(s$Duration__in_seconds_, 0.2)),
      consumer = sum(consumer_vig),
      strictconseq = sum(wtp_ratio_class == 0, na.rm=T),
      weakconseq = sum(wtp_ratio < 1, na.rm=T)
    ) %>%
    mutate(
      freq_main = n / sum(n, na.rm=T),
      freq_wgt = wgt / sum(wgt, na.rm=T),
      freq_nospeeders = nospeeder / sum(nospeeder),
      freq_consumer = consumer / sum(consumer),
      freq_strictconseq = strictconseq / sum(strictconseq),
      freq_weakconseq = weakconseq / sum(weakconseq)
    ) %>% filter(code_effect_simple %in% c("dampening", "tiny")),
  s %>% group_by(code_effect_simple) %>%
    summarise(
      n = n(),
      wgt = sum(wgt),
      nospeeder = sum(Duration__in_seconds_ > quantile(s$Duration__in_seconds_, 0.2)),
      consumer = sum(consumer_vig),
      strictconseq = sum(wtp_ratio_class == 0, na.rm=T),
      weakconseq = sum(wtp_ratio < 1, na.rm=T)
    ) %>%
    mutate(
      vig = "pooled",
      freq_main = n / sum(n, na.rm=T),
      freq_wgt = wgt / sum(wgt, na.rm=T),
      freq_nospeeders = nospeeder / sum(nospeeder),
      freq_consumer = consumer / sum(consumer),
      freq_strictconseq = strictconseq / sum(strictconseq),
      freq_weakconseq = weakconseq / sum(weakconseq)
    ) %>% filter(code_effect_simple != "dampening"),
  s %>% group_by(code_effect_simple=code_effect) %>%
    summarise(
      n = n(),
      wgt = sum(wgt),
      nospeeder = sum(Duration__in_seconds_ > quantile(s$Duration__in_seconds_, 0.2)),
      consumer = sum(consumer_vig),
      strictconseq = sum(wtp_ratio_class == 0, na.rm=T),
      weakconseq = sum(wtp_ratio < 1, na.rm=T)
    ) %>%
    mutate(
      vig = "pooled",
      freq_main = n / sum(n, na.rm=T),
      freq_wgt = wgt / sum(wgt, na.rm=T),
      freq_nospeeders = nospeeder / sum(nospeeder),
      freq_consumer = consumer / sum(consumer),
      freq_strictconseq = strictconseq / sum(strictconseq),
      freq_weakconseq = weakconseq / sum(weakconseq)
    ) %>% filter(code_effect_simple %in% c("dampening", "tiny"))
)
df = df %>% pivot_longer(starts_with("freq_"), values_to="freq")
df$name = gsub("^freq_", "", df$name)
df$code_effect_simple = factor(df$code_effect_simple, levels=c("other", "spillover", "naive", "tiny", "dampening"))

# Don't display other fraction
df = df %>% filter(code_effect_simple != "other")

# Order
vigs = c(
  "pooled" = "All cases",
  "fuel" = "Fuel",
  "meat" = "Meat",
  "energy" = "Energy",
  "flights" = "Plane trips",
  "subelectricity" = "Green\nelectricity",
  "subhouse" = "Energy-eff.\nhome",
  "subclothing" = "2nd-hand\nclothing",
  "subcoffee" = "Fairtrade\ncoffee"
)
df$vig = factor(df$vig, levels=names(vigs))

# Robustness tests
tests = c(
  "main" = 1,
  "wgt" = 2,
  "nospeeders" = 3,
  "consumer" = 4,
  "strictconseq" = 5,
  "weakconseq" = 6
)


######################
#### PLOT ############
######################
plots = list()

plots$fig_explanations_robust = ggplot(df, aes(x = name, y = freq, fill = code_effect_simple)) +
  geom_bar(stat = "identity") +
  facet_grid(.~vig, labeller = as_labeller(vigs)) +
  xlab("Specification") + ylab("Frequency") +
  coord_cartesian(ylim=c(0, 1)) +
  scale_x_discrete(limits = names(tests), labels=tests) +
  scale_y_continuous(breaks=seq(0, 1, 0.2), labels=paste0(100*seq(0, 1, 0.2), "%")) +
  scale_fill_manual(
    values=c("darkolivegreen2", "skyblue", "red1", "red3"),
    labels=c("Positive multiplier  ", "One-to-one impact  ", "Dampened impact:   \nToo small to matter", "Dampened impact:   \nOffset by others   "),
    name = "**Explanation**: Respondents explain a ..."
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "rob_con"),
    plot.title = element_text(size=18, face="bold", hjust=0.5),
    legend.position = "bottom",
    legend.direction = "vertical",
    legend.title = element_markdown(size=14, hjust=0.5),
    #legend.spacing.y = unit(0, 'cm'),
    panel.grid.major.x = element_blank(),
    panel.grid.minor = element_blank(),
    strip.text = element_text(size=14, face="bold"),
    axis.title = element_text(size=14, face="bold"),
    axis.text.x = element_markdown(size=14),
    axis.text.y = element_text(size=14),
    legend.text = element_text(size=14)
  ) +
  guides(fill = guide_legend(nrow = 1, reverse = TRUE))

# Save generated graph.
pdf(path_to("figures", "fig_explanations_robust.pdf"), width=11, height=5.5)
print(plots$fig_explanations_robust)
dev.off()  

